using Gadfly
using DataFrames
using Colors
# E.g.
set_default_plot_size(26cm, 14cm)
theta1(x) = float(x)^2*(1.0-x)^2
theta2(x) = (x >= 0.25 && x <= 0.75) ? 1 : 0
intervall_discretisation(n) = linspace(1.0/n, 1 - 1.0/n, n-1);
function gradient_descent(f0, A, Theta; gamma=0.0, eps=1.0e-6)
k = 1
f = f0
grad = A\(A\f-Theta) + gamma * f
while (norm(grad) > eps) && (k < 2e4)
k = k + 1
grad = A\(A\f-Theta) + gamma * f
A_inv_grad = A\grad
grad_dot = dot(grad, grad)
sigma_exact = grad_dot / (dot(A_inv_grad, A_inv_grad) + gamma * grad_dot)
f = f - sigma_exact * grad
end
return f, k
end
function waerme(theta, n; eps=1.0e-6, gamma=0.0)
e = ones(n-1)
A = Tridiagonal(-e[2:end], 2e, -e[2:end]) * n^2
Theta = map(theta, intervall_discretisation(n))
f, k = gradient_descent(e, A, Theta, gamma=gamma, eps=eps)
b = A\f
return Dict("f" => f, "b" => b, "k" => k, "Theta" => Theta)
end
Berechne für verschiedene $\gamma$ und $n$
#gamma_values = [0; [(iseven(x) ? 1.0 : 0.5)*(10.0^(-floor(x / 2))) for x=16:-1:8]];
gamma_values = [0; 5.0e-8; 1.0e-7; 1.0e-6; 1.0e-5; 1.0e-4; 1.0e-3];
print(gamma_values)
m = length(gamma_values) - 1
cmap = colormap("Blues", m + 1)[end:-1:2]
values_theta1_10 = [waerme(theta1, 10, gamma=x) for x in gamma_values];
values_theta1_100 = [waerme(theta1, 100, gamma=x) for x in gamma_values];
values_theta1_1000 = [waerme(theta1, 1000, gamma=x) for x in gamma_values];
values_theta2_10 = [waerme(theta2, 10, gamma=x) for x in gamma_values];
values_theta2_100 = [waerme(theta2, 100, gamma=x) for x in gamma_values];
values_theta2_1000 = [waerme(theta2, 1000, gamma=x) for x in gamma_values];
set_default_plot_size(25.6cm, 19.6cm)
function plot_control_gammas(gamma_values ;theta="theta1", n=1000)
values
if (theta == "theta1")
values = n == 1000 ? values_theta1_1000 : ((n == 100) ? values_theta1_100 : values_theta1_10)
else
values = n == 1000 ? values_theta2_1000 : ((n == 100) ? values_theta2_100 : values_theta2_10)
end
m = length(gamma_values) - 1
cmap = colormap("Blues", m + 1)[end:-1:2]
legendValues = [
"Without regularization & t = "string(values[1]["k"]);
map(x -> "𝛾 = "string(x[1])" & t = "string(x[2]["k"]), zip(gamma_values[2:end], values[2:end]))
]
legendColors = [colorant"red"; cmap]
layers = [layer(
x=intervall_discretisation(n),
y=values[i]["f"],
Geom.line,
Theme(default_color=legendColors[i], line_width=0.5mm)
) for i=1:m]
return plot(layers...,
Guide.manual_color_key("Legend", legendValues, legendColors),
Guide.title("Optimal control for " * (theta == "theta1" ? "θ₁" : "θ₂") * " and N=" * string(n)),
Guide.YLabel("Optimal control fN")
)
end
function plot_solution_gammas(gamma_values ;theta="theta1", n=1000)
values
if (theta == "theta1")
values = n == 1000 ? values_theta1_1000 : ((n == 100) ? values_theta1_100 : values_theta1_10)
else
values = n == 1000 ? values_theta2_1000 : ((n == 100) ? values_theta2_100 : values_theta2_10)
end
m = length(gamma_values) - 1
cmap = colormap("Blues", m + 1)[end:-1:2]
legendValues = [
"Required temperatur";
"Without regularization & t = "string(values[1]["k"]);
map(x -> "𝛾 = "string(x[1])" & t = "string(x[2]["k"]), zip(gamma_values[2:end], values[2:end]))
]
legendColors = [colorant"green"; colorant"red"; cmap]
layers =
[layer(
x = intervall_discretisation(n),
y = values[1]["Theta"],
Geom.point,
Theme(default_color=colorant"green", highlight_width=0.05mm, default_point_size=0.65mm)
);
[layer(
x=intervall_discretisation(n),
y=values[i]["b"],
Geom.line,
Theme(default_color=legendColors[i+1], line_width=0.5mm)
) for i=1:m]]
return plot(layers...,
Guide.manual_color_key("Legend", legendValues, legendColors),
Guide.title("Optimal solution yN; for " * (theta == "theta1" ? "θ₁" : "θ₂") * " and N=" * string(n)),
Guide.YLabel("State yN")
)
end
hstack(
vstack(
plot_control_gammas(gamma_values, theta="theta1", n=10),
plot_control_gammas(gamma_values, theta="theta1", n=100),
plot_control_gammas(gamma_values, theta="theta1", n=1000)
),
vstack(
plot_solution_gammas(gamma_values, theta="theta1", n=10),
plot_solution_gammas(gamma_values, theta="theta1", n=100),
plot_solution_gammas(gamma_values, theta="theta1", n=1000)
))
hstack(
vstack(
plot_control_gammas(gamma_values, theta="theta2", n=10),
plot_control_gammas(gamma_values, theta="theta2", n=100),
plot_control_gammas(gamma_values, theta="theta2", n=1000)
),
vstack(
plot_solution_gammas(gamma_values, theta="theta2", n=10),
plot_solution_gammas(gamma_values, theta="theta2", n=100),
plot_solution_gammas(gamma_values, theta="theta2", n=1000)
))
vstack(
plot_control_gammas(gamma_values, theta="theta1", n=10),
plot_control_gammas(gamma_values, theta="theta1", n=100),
plot_control_gammas(gamma_values, theta="theta1", n=1000)
)
vstack(
plot_solution_gammas(gamma_values, theta="theta1", n=10),
plot_solution_gammas(gamma_values, theta="theta1", n=100),
plot_solution_gammas(gamma_values, theta="theta1", n=1000)
)
vstack(
plot_control_gammas(gamma_values, theta="theta2", n=10),
plot_control_gammas(gamma_values, theta="theta2", n=100),
plot_control_gammas(gamma_values, theta="theta2", n=1000)
)
vstack(
plot_solution_gammas(gamma_values, theta="theta2", n=10),
plot_solution_gammas(gamma_values, theta="theta2", n=100),
plot_solution_gammas(gamma_values, theta="theta2", n=1000)
)
vstack(
plot_control_gammas(gamma_values, theta="theta1", n=100),
plot_solution_gammas(gamma_values, theta="theta1", n=100),
)
vstack(
plot_control_gammas(gamma_values, theta="theta2", n=100),
plot_solution_gammas(gamma_values, theta="theta2", n=100),
)